自己相関行列の逆の観察

自己相関行列の逆行列の観察

なんといってもまずは観察から。

In [1]:
import math
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

NUM_SAMPLES = 4096

# データ生成
sindata = np.zeros(NUM_SAMPLES)
gaussnoise = np.random.normal(0, 1, NUM_SAMPLES)
laplacenoise = np.random.laplace(0, 1, NUM_SAMPLES)
for i in range(NUM_SAMPLES):
    sindata[i] = math.sin(2 * math.pi *  440 * i / 48000.0)
In [2]:
# 自己相関行列を計算
def calc_auto_corration_matrix(data, dim):
    S = np.zeros((dim, dim))
    for smpl in range(dim, data.size, 1):
        xvec = data[smpl - dim:smpl].reshape((dim, 1))
        S += np.dot(xvec, xvec.T)
    return S / (data.size - dim)

# 分散1, 平均0化
def normalize(data):
    return (data - np.mean(data)) / np.std(data)

# 自己相関行列とその逆の表示
def plot_gt_auto_correlation(data, dim):
    S = calc_auto_corration_matrix(data, dim)
    _, axs = plt.subplots(ncols=2, figsize=(10, 4))
    sns.heatmap(S, ax=axs[0])
    axs[0].set_title('Auto Correlation (Grand Trues)')
    sns.heatmap(np.linalg.inv(S), ax=axs[1])
    axs[1].set_title('Inverse Auto Correlation (Grand Trues)')
In [3]:
# ノイズの比をいじって調べる
plot_gt_auto_correlation(normalize(gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 1.0  * gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.75 * gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.5  * gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.25 * gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.0  * gaussnoise), 16)

考察

  • ノイズが乗って初めて対角要素が強調される。ノイズが乗ることが大事。
    • ある種当然か。ノイズがあるから揺らぎが生まれて分散が出る。
  • 自己相関行列の逆、ノイズがない場合を除いて、定数倍の違いがあるだけでは？色のパターンがほぼ同じ。

次に調べたいこと

  • ノイズにも相関がある場合は？
In [4]:
# 信号に1次の相関を付与
def make_1st_correlation(data, corr):
    corr_data = data.copy()
    for smpl in range(1, data.size, 1):
        corr_data[smpl] += corr * corr_data[smpl - 1]
    return corr_data

corr_gaussnoise = make_1st_correlation(gaussnoise, 0.95)
plot_gt_auto_correlation(normalize(corr_gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 1.0  * corr_gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.75 * corr_gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.5  * corr_gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.25 * corr_gaussnoise), 16)
plot_gt_auto_correlation(normalize(sindata + 0.0  * corr_gaussnoise), 16)

相関付きノイズ・考察

  • 3重対角要素が発生した。
    • 対角にある0次相関と比べると、3重対角要素の1次相関は負に振り切った値
    • それ意外の要素はほぼ0では…？
  • 1行目と最後の行を除いて巡回行列になっていないか？
    • いや、ちょっと正確じゃない。行列を対角線で切り分けた時、上三角に置かれている要素を展開しているだけに見える

次に見たいこと

  • SN比との関係
    • いよいよ正弦波入力が関係ないように見えてきた。ノイズだけが本質になりそうな予感がする。
    • 定数倍を除いて皆同じになるか？
In [5]:
for ratio in np.linspace(0.0, 1.0, 10):
    mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(gaussnoise)
    mixed = normalize(mixed)
    plot_gt_auto_correlation(mixed, 16)
    
for ratio in np.linspace(0.0, 1.0, 10):
    mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(corr_gaussnoise)
    mixed = normalize(mixed)
    plot_gt_auto_correlation(mixed, 16)
In [6]:
# 自己相関行列の逆行列の比較: (data+noise, noiseのみ)の比較
def plot_compare_inverse_auto_correlation(mixed, noise, dim):
    Smix   = calc_auto_corration_matrix(mixed, dim)
    Snoise = calc_auto_corration_matrix(noise, dim)
    Smix_inv = np.linalg.inv(Smix)
    Snoise_inv = np.linalg.inv(Snoise)
    _, axs = plt.subplots(ncols=3, figsize=(14, 4))
    sns.heatmap(Smix_inv, ax=axs[0])
    axs[0].set_title('Inverse Auto Correlation for Data+Noise')
    sns.heatmap(Snoise_inv, ax=axs[1])
    axs[1].set_title('Inverse Auto Correlation for Noise')
    sns.heatmap(np.abs(Smix_inv - Snoise_inv), ax=axs[2])
    axs[2].set_title('Absolute Error')
    
for ratio in np.linspace(0.0, 1.0, 10):
    mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(gaussnoise)
    mixed = normalize(mixed)
    plot_compare_inverse_auto_correlation(mixed, normalize(gaussnoise), 16)
    
for ratio in np.linspace(0.0, 1.0, 10):
    mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(corr_gaussnoise)
    mixed = normalize(mixed)
    plot_compare_inverse_auto_correlation(mixed, normalize(corr_gaussnoise), 16)
In [7]:
# 自己相関行列の逆行列の比較: (data+noise, noiseのみ)の比較 L2ノルム正規化して比較
def plot_norm_compare_inverse_auto_correlation(mixed, noise, dim):
    Smix   = calc_auto_corration_matrix(mixed, dim)
    Snoise = calc_auto_corration_matrix(noise, dim)
    Smix_inv = np.linalg.inv(Smix)
    Smix_inv /= np.linalg.norm(Smix_inv)
    Snoise_inv = np.linalg.inv(Snoise)
    Snoise_inv /= np.linalg.norm(Snoise_inv)
    _, axs = plt.subplots(ncols=3, figsize=(14, 4))
    sns.heatmap(Smix_inv, ax=axs[0])
    axs[0].set_title('Inverse Auto Correlation for Data+Noise')
    sns.heatmap(Snoise_inv, ax=axs[1])
    axs[1].set_title('Inverse Auto Correlation for Noise')
    sns.heatmap(20 * np.log10(np.abs(Smix_inv - Snoise_inv)), ax=axs[2])
    axs[2].set_title('Absolute Error [dB]')
    
for ratio in np.linspace(0.0, 1.0, 10):
    mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(gaussnoise)
    mixed = normalize(mixed)
    plot_norm_compare_inverse_auto_correlation(mixed, normalize(gaussnoise), 16)
    
for ratio in np.linspace(0.0, 1.0, 10):
    mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(corr_gaussnoise)
    mixed = normalize(mixed)
    plot_norm_compare_inverse_auto_correlation(mixed, normalize(corr_gaussnoise), 16)

<ipython-input-7-5edd95170ad1>:14: RuntimeWarning: divide by zero encountered in log10
  sns.heatmap(20 * np.log10(np.abs(Smix_inv - Snoise_inv)), ax=axs[2])

考察

  • ノルムを揃えたら、正弦波+ノイズとノイズがほぼ一致している。
    • 入力データは定数倍にしか効いてこない？
    • 相関付きノイズの場合はより誤差が小さい。
  • 相関付きノイズは入力の大きさに関わらず、逆行列が一定に見える。
    • 相関なしのi.i.d.ノイズでは入力音声のパワーに応じて対角要素以外が明るくなる。
  • n重対角要素の誤差が大きい
    • n重対角要素外では誤差が殆どない（暗い）。定数倍除いてもほぼ誤差なし。

次に気になること

  • n重対角要素以外、ほぼ0では？
    • 要素の絶対値をプロットしてみる
    • これが分かると、スパース性を使える。
In [8]:
# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_abs_inverse_auto_correlation(data, dim):
    S = calc_auto_corration_matrix(data, dim)
    Sinv = np.linalg.inv(S)
    Sinv /= np.linalg.norm(Sinv)
    _, axs = plt.subplots(ncols=2, figsize=(10, 4))
    sns.heatmap(np.abs(Sinv), ax=axs[0])
    axs[0].set_title('Absolute Inverse Auto Correlation')
    sns.heatmap(20 * np.log10(np.abs(Sinv)), ax=axs[1])
    axs[1].set_title('Log-Absolute Inverse Auto Correlation [dB]')
    
for ratio in np.linspace(0.0, 1.0, 10):
    mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(gaussnoise)
    mixed = normalize(mixed)
    plot_abs_inverse_auto_correlation(mixed, 16)
    
for ratio in np.linspace(0.0, 1.0, 10):
    mixed = ratio * normalize(sindata) + (1.0 - ratio) * normalize(corr_gaussnoise)
    mixed = normalize(mixed)
    plot_abs_inverse_auto_correlation(mixed, 16)

考察

  • n重対角要素を除き、0と見て良いと思える。
    • 少なくとも近似的には言えるはず。相関付きノイズではさらにはっきりと言える。
    • 0とみなした要素は全て計算を省いてよさそう。
    • さっきも触れたけど、行列に対角線を引いた時の上三角要素を全てに展開しているように見える。

次にやること

  • 近似計算の影響は？
    • 最終的には途中経過が見たいけど、まずは最後にどんな結果が得られているか見よう。
In [9]:
# Sherman–Morrison法による自己相関の逆の逐次計算
def calc_inv_auto_corr_by_sm(data, dim, forget):
    Sinv = np.eye(dim)
    for smpl in range(dim, data.size, 1):
        xvec = data[smpl - dim:smpl].reshape((dim, 1))
        Sx = np.dot(Sinv, xvec)
        Sxxt = np.dot(Sx, Sx.T) 
        Sinv -= Sxxt / (forget + np.dot(xvec.T, Sx))
        Sinv /= forget
    return Sinv

# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_sm_and_gt(data, dim):
    S = calc_auto_corration_matrix(data, dim)
    Sinv = np.linalg.inv(S)
    Sinv_sm = calc_inv_auto_corr_by_sm(data, dim, 0.98)
    Sinv /= np.linalg.norm(Sinv)
    Sinv_sm /= np.linalg.norm(Sinv_sm)
    _, axs = plt.subplots(ncols=3, figsize=(15, 4))
    sns.heatmap(Sinv, ax=axs[0])
    axs[0].set_title('Inverse Auto Correlation(Grand Trues)')
    sns.heatmap(Sinv_sm, ax=axs[1])
    axs[1].set_title('Inverse Auto Correlation(Sherman-Morrison)')
    sns.heatmap(np.abs(Sinv - Sinv_sm), ax=axs[2])
    axs[2].set_title('Absolute Error')
    
plot_auto_corr_compare_sm_and_gt(gaussnoise, 16)
plot_auto_corr_compare_sm_and_gt(normalize(corr_gaussnoise), 16)

# ノイズ以外の要素は、もはや気にしない！（上の考察より）
# plot_auto_corr_compare_sm_and_gt(normalize(sindata + gaussnoise), 16)
# plot_auto_corr_compare_sm_and_gt(normalize(sindata + corr_gaussnoise), 16)

考察

  • (n重)非対角要素の誤差が大きい。斜め方向に誤差が出ている。
    • (n重)非対角成分は0に落としちゃっていい。スパース制約が大事に見える。

次に見たいもの

  • Graphical LASSOとの比較。
In [10]:
from sklearn.covariance import GraphicalLasso
import scipy as sp

# Graphical Lassoによる精度行列計算（平均0だから自己相関行列の逆）
def calc_inv_auto_corr_by_gl(data, dim, alpha):
    X = np.zeros((data.size - dim, dim))
    for smpl in range(dim, data.size, 1):
        X[smpl - dim:] = data[smpl - dim:smpl]
    gl = GraphicalLasso(alpha=alpha)
    gl.fit(X)
    return gl.precision_

# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_gl_and_gt(data, dim, alpha):
    S = calc_auto_corration_matrix(data, dim)
    Sinv = np.linalg.inv(S)
    Sinv_gl = calc_inv_auto_corr_by_gl(data, dim, alpha)
    Sinv /= np.linalg.norm(Sinv)
    Sinv_gl /= np.linalg.norm(Sinv_gl)
    _, axs = plt.subplots(ncols=3, figsize=(15, 4))
    sns.heatmap(Sinv, ax=axs[0])
    axs[0].set_title('Inverse Auto Correlation(Grand Trues)')
    sns.heatmap(Sinv_gl, ax=axs[1])
    axs[1].set_title('Inverse Auto Correlation(Graphical LASSO alpha: ' + str(alpha) + ')')
    sns.heatmap(np.abs(Sinv - Sinv_gl), ax=axs[2])
    axs[2].set_title('Absolute Error')

plot_auto_corr_compare_gl_and_gt(gaussnoise, 16, 0.1)
plot_auto_corr_compare_gl_and_gt(corr_gaussnoise, 16, 0.1)

考察

  • 残差を見ると、Sherman-Morrisonよりも良さそうに見える。
  • いよいよ、実装してみるべきなのかもしれない。
  • 実データに対しても同様の傾向なのかは見ておきたい。

成果報告向けプロット

In [11]:
# 真のRの逆行列はどうなっているか？
plot_gt_auto_correlation(normalize(gaussnoise), 16)
plt.savefig("r_invr_gt_whitenoise.png")
plot_gt_auto_correlation(normalize(corr_gaussnoise), 16)
plt.savefig("r_invr_gt_corrnoise.png")
# Sherman-Morrisonとの比較
plot_auto_corr_compare_sm_and_gt(gaussnoise, 16)
plt.savefig("invr_comp_sm_and_gt_whitenoise.png")
plot_auto_corr_compare_sm_and_gt(normalize(corr_gaussnoise), 16)
plt.savefig("invr_comp_sm_and_gt_corrnoise.png")
# glassoとの比較
plot_auto_corr_compare_gl_and_gt(gaussnoise, 16, 0.1)
plt.savefig("invr_comp_gl_and_gt_whitenoise.png")
plot_auto_corr_compare_gl_and_gt(corr_gaussnoise, 16, 0.1)
plt.savefig("invr_comp_gl_and_gt_corrnoise.png")

追試：RLSの更新則で見ている自己相関行列

  • Sharman-Morrisonで逐次計算しないやり方に該当する。観察。
  • 性能としてはRLSの方が良いことが分かっている。
  • 上とは何か違うものを見ている可能性が大。
In [12]:
# RLSが考えている自己相関行列の計算
def calc_inv_auto_corr_by_rls(data, dim, forget):
    S = np.eye(dim)
    for smpl in range(dim, data.size, 1):
        xvec = data[smpl - dim:smpl].reshape((dim, 1))
        S = forget * S + xvec @ xvec.T
        # S = (smpl * S + xvec @ xvec.T) / (smpl + 1)
    return S, np.linalg.inv(S)

# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_rls_and_gt(data, dim, forget):
    S = calc_auto_corration_matrix(data, dim)
    Sinv = np.linalg.inv(S)
    Srls, Sinv_rls = calc_inv_auto_corr_by_rls(data, dim, forget)
    Sinv /= np.linalg.norm(Sinv)
    Sinv_rls /= np.linalg.norm(Sinv_rls)
    _, axs = plt.subplots(ncols=3, figsize=(15, 4))
    # sns.heatmap(Sinv, ax=axs[0])
    # axs[0].set_title('Inverse Auto Correlation(Grand Trues)')
    sns.heatmap(Srls, ax=axs[0])
    axs[0].set_title('Auto Correlation(RLS) Num Samples:%d' % data.size)
    sns.heatmap(Sinv_rls, ax=axs[1])
    axs[1].set_title('Inverse Auto Correlation(RLS)')
    sns.heatmap(np.abs(Sinv - Sinv_rls), ax=axs[2])
    axs[2].set_title('Absolute Error between the GT')
    
gaussnoise = np.random.normal(0, 1, NUM_SAMPLES)
corr_gaussnoise = make_1st_correlation(gaussnoise, 0.95)

for nsmpls in np.arange(100, 1000, 100):
    plot_auto_corr_compare_rls_and_gt(gaussnoise[0:nsmpls], 16, 0.999)
    plot_auto_corr_compare_rls_and_gt(corr_gaussnoise[0:nsmpls], 16, 0.999)

所感

  • ほぼリファレンス(GT)と変わらないように見えるのだが...
    • 各サンプルでGTに近ければ良い結果は出るはず。
  • 定常状態(GTで得られる答え)との誤差は十分に小さいように見受けられる。
  • 自分が考えてるAR(1)仮定をぶち込んで見てみよう。
In [13]:
# AR(1)仮定による逆行列計算
def calc_inv_auto_corr_by_ar1(data, dim, forget):
    var00 = 1.0
    var01 = 0.0
    tri_mask = np.multiply(np.tri(dim, dim, -1), np.tri(dim, dim, 1).T) 
    tri_mask += tri_mask.T
    # for smpl in range(dim, data.size, 1):
        # var00 = forget * var00 + (1.0 - forget) * data[smpl] ** 2.0
        # var01 = forget * var01 + (1.0 - forget) * data[smpl] * data[smpl - 1]
    var00 = np.sum(data ** 2) / data.size
    var01 = np.sum(data[0:data.size - 1] * data[1:data.size]) / (data.size - 1)
    cor01 = var01 / var00
    Rinv = (1.0 + (cor01 ** 2)) * np.eye(dim)
    Rinv[0, 0] = Rinv[dim-1, dim-1] = 1.0
    Rinv += -cor01 * tri_mask
    Rinv /= var00 * (1.0 - (cor01 ** 2))
    return Rinv

# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_ar1_and_gt(data, dim, forget):
    S = calc_auto_corration_matrix(data, dim)
    Sinv = np.linalg.inv(S)
    Sinv_ar1 = calc_inv_auto_corr_by_ar1(data, dim, forget)
    Sinv_gl = calc_inv_auto_corr_by_gl(data, dim, 0.2)
    Sinv /= np.linalg.norm(Sinv)
    Sinv_ar1 /= np.linalg.norm(Sinv_ar1)
    Sinv_gl /= np.linalg.norm(Sinv_gl)
    _, axs = plt.subplots(ncols=3, figsize=(15, 4))
    sns.heatmap(Sinv, ax=axs[0])
    axs[0].set_title('Inverse Auto Correlation(Grand Trues)')
    # sns.heatmap(Sinv_gl, ax=axs[0])
    # axs[0].set_title('Inverse Auto Correlation(glasso)')
    sns.heatmap(Sinv_ar1, ax=axs[1])
    axs[1].set_title('Inverse Auto Correlation(AR(1))')
    sns.heatmap(np.abs(Sinv - Sinv_ar1), ax=axs[2])
    axs[2].set_title('Absolute Error')

for nsmpls in np.arange(30, 200, 20):
    plot_auto_corr_compare_ar1_and_gt(gaussnoise[0:nsmpls], 16, 0.999)
    plot_auto_corr_compare_ar1_and_gt(corr_gaussnoise[0:nsmpls], 16, 0.999)
In [29]:
# AR(2)仮定による逆行列計算
def calc_inv_auto_corr_by_ar2(data, dim, forget):
    var00 = 1.0
    var01 = 0.0
    var02 = 0.0
    tri_mask_ord1 = np.multiply(np.tri(dim, dim, -1), np.tri(dim, dim, 1).T) 
    tri_mask_ord1 += tri_mask_ord1.T
    tri_mask_ord2 = np.multiply(np.tri(dim, dim, -2), np.tri(dim, dim, 2).T) 
    tri_mask_ord2 += tri_mask_ord2.T
    # for smpl in range(dim, data.size, 1):
        # var00 = forget * var00 + (1.0 - forget) * data[smpl] ** 2.0
        # var01 = forget * var01 + (1.0 - forget) * data[smpl] * data[smpl - 1]
    var00 = np.sum(data ** 2) / data.size
    var01 = np.sum(data[0:data.size - 1] * data[1:data.size]) / (data.size - 1)
    var02 = np.sum(data[0:data.size - 2] * data[2:data.size]) / (data.size - 2)
    cor01 = var01 / var00
    cor02 = var02 / var00
    phi1 = cor01 * (1.0 - cor02) / (1.0 - (cor01 ** 2))
    phi2 = (cor02 - (cor01 ** 2)) / (1.0 - (cor01 ** 2))
    diag = (1.0 + (phi1 ** 2) + (phi2 ** 2)) * np.eye(dim)
    diag[0, 0] = diag[dim-1, dim-1] = 1.0
    diag[1, 1] = diag[dim-2, dim-2] = 1.0 + (phi1 ** 2)
    tri_mask_ord1 = -phi1 * (1.0 - phi2) * tri_mask_ord1
    tri_mask_ord1[0, 1] = tri_mask_ord1[1, 0] \
        = tri_mask_ord1[dim-1, dim-2] = tri_mask_ord1[dim-2, dim-1] = -phi1
    tri_mask_ord2 = -phi2 * tri_mask_ord2
    Rinv = diag + tri_mask_ord1 + tri_mask_ord2
    Rinv /= var00 * (1.0 + phi2) * (((1.0 - phi2) ** 2) - (phi1 ** 2)) / (1.0 - phi2)
    return Rinv

# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_ar2_and_gt(data, dim, forget):
    S = calc_auto_corration_matrix(data, dim)
    Sinv = np.linalg.inv(S)
    Sinv_ar2 = calc_inv_auto_corr_by_ar2(data, dim, forget)
    Sinv_gl = calc_inv_auto_corr_by_gl(data, dim, 0.1)
    Sinv /= np.linalg.norm(Sinv)
    Sinv_ar2 /= np.linalg.norm(Sinv_ar2)
    Sinv_gl /= np.linalg.norm(Sinv_gl)
    _, axs = plt.subplots(ncols=3, figsize=(15, 4))
    sns.heatmap(Sinv, ax=axs[0])
    axs[0].set_title('Inverse Auto Correlation(Grand Trues)')
    # sns.heatmap(Sinv_gl, ax=axs[0])
    # axs[0].set_title('Inverse Auto Correlation(glasso)')
    sns.heatmap(Sinv_ar2, ax=axs[1])
    axs[1].set_title('Inverse Auto Correlation(AR(2))')
    sns.heatmap(np.abs(Sinv - Sinv_ar2), ax=axs[2])
    axs[2].set_title('Absolute Error between the glasso')

# for nsmpls in np.arange(30, 200, 20):
for nsmpls in np.arange(100, 1000, 200):
    plot_auto_corr_compare_ar2_and_gt(gaussnoise[0:nsmpls], 16, 0.999)
    plot_auto_corr_compare_ar2_and_gt(corr_gaussnoise[0:nsmpls], 16, 0.999)

/usr/local/lib/python3.8/site-packages/sklearn/covariance/_graph_lasso.py:264: ConvergenceWarning: graphical_lasso: did not converge after 100 iteration: dual gap: -1.048e-02
  warnings.warn('graphical_lasso: did not converge after '
/usr/local/lib/python3.8/site-packages/sklearn/covariance/_graph_lasso.py:264: ConvergenceWarning: graphical_lasso: did not converge after 100 iteration: dual gap: 6.096e-03
  warnings.warn('graphical_lasso: did not converge after '

所感

  • AR(1), AR(2)（今実装した）は誤差がそれなりに(0.02くらい)出ている。
    • RLSの考え方、つまりSherman-Morrisonの誤差より大きい。
  • スパースなglassoで求めた結果と比較しても傾向は同じ。
In [15]:
# 自己相関行列の逆の絶対値、対数絶対値の表示
def plot_auto_corr_compare_ar1_ar2_gt(data, dim, forget):
    S = calc_auto_corration_matrix(data, dim)
    Sinv = np.linalg.inv(S)
    Sinv_ar1 = calc_inv_auto_corr_by_ar1(data, dim, forget)
    Sinv_ar2 = calc_inv_auto_corr_by_ar2(data, dim, forget)
    _, Sinv_rls = calc_inv_auto_corr_by_rls(data, dim, forget)
    Sinv /= np.linalg.norm(Sinv)
    Sinv_ar1 /= np.linalg.norm(Sinv_ar1)
    Sinv_ar2 /= np.linalg.norm(Sinv_ar2)
    Sinv_rls /= np.linalg.norm(Sinv_rls)
    _, axs = plt.subplots(ncols=5, figsize=(16, 3))
    sns.heatmap(Sinv_ar1, ax=axs[0])
    axs[0].set_title('Inverse Auto Correlation(AR(1))')
    sns.heatmap(Sinv_ar2, ax=axs[1])
    axs[1].set_title('Inverse Auto Correlation(AR(2))')
    sns.heatmap(np.abs(Sinv - Sinv_ar1), ax=axs[2], vmin=0, vmax=0.1)
    axs[2].set_title('Absolute Error (AR(1))')
    sns.heatmap(np.abs(Sinv - Sinv_ar2), ax=axs[3], vmin=0, vmax=0.1)
    axs[3].set_title('Absolute Error (AR(2))')
    sns.heatmap(np.abs(Sinv - Sinv_rls), ax=axs[4], vmin=0, vmax=0.1)
    axs[4].set_title('Absolute Error (RLS)')
In [16]:
# i.i.d. に対して実験
for nsmpls in np.arange(50, 600, 50):
    print(nsmpls)
    plot_auto_corr_compare_ar1_ar2_gt(gaussnoise[0:nsmpls], 16, 0.999)
    plt.show()

50

100

150

200

250

300

350

400

450

500

550
In [17]:
# 相関付きに対して実験
for nsmpls in np.arange(50, 600, 50):
    print(nsmpls)
    plot_auto_corr_compare_ar1_ar2_gt(corr_gaussnoise[0:nsmpls], 16, 0.999)
    plt.show()

50

100

150

200

250

300

350

400

450

500

550

所感

  • RLSはかなり正確な予測をしており、GTへの収束も速い。
  • 500サンプルくらいまで行けばほぼGTと一致する。が、その間に勝負が決まってしまうっぽい。
  • 入力が最大で1サンプル相関だからというのがあるが、AR(1)とAR(2)であまり性能がかわらない。

次は…

  • GTとのRMSE曲線書いてみると良さそう→今すぐやってみる
  • 共役勾配法をRLS流でやればうまくいくはず。自己相関行列を逐次正しく求めているから。早くやるのだ。
In [18]:
# 真の逆行列との誤差RMSEの変化をプロット
def plot_inv_auto_corr_rmse_curve(dim, forget, smpl_range, data_generator, num_trials):
    num_ranges = smpl_range.size
    list_ar1 = np.zeros(num_ranges)
    list_ar2 = np.zeros(num_ranges)
    list_rls = np.zeros(num_ranges)
    list_gl = np.zeros(num_ranges)
    for _ in range(num_trials):
        for idx, nsmpl in enumerate(smpl_range):
            data = data_generator(nsmpl)
            S = calc_auto_corration_matrix(data, dim)
            Sinv = np.linalg.inv(S)
            Sinv_ar1 = calc_inv_auto_corr_by_ar1(data, dim, forget)
            Sinv_ar2 = calc_inv_auto_corr_by_ar2(data, dim, forget)
            _, Sinv_rls = calc_inv_auto_corr_by_rls(data, dim, forget)
            Sinv_gl = calc_inv_auto_corr_by_gl(data, dim, 1e-1)
            Sinv /= np.linalg.norm(Sinv)
            Sinv_ar1 /= np.linalg.norm(Sinv_ar1)
            Sinv_ar2 /= np.linalg.norm(Sinv_ar2)
            Sinv_rls /= np.linalg.norm(Sinv_rls)
            Sinv_gl /= np.linalg.norm(Sinv_gl)
            list_ar1[idx] += np.linalg.norm(Sinv - Sinv_ar1)
            list_ar2[idx] += np.linalg.norm(Sinv - Sinv_ar2)
            list_rls[idx] += np.linalg.norm(Sinv - Sinv_rls)
            list_gl[idx] += np.linalg.norm(Sinv - Sinv_gl)
    list_ar1 /= num_trials
    list_ar2 /= num_trials
    list_rls /= num_trials
    list_gl /= num_trials
    plt.plot(smpl_range, list_ar1, label='AR(1)')
    plt.plot(smpl_range, list_ar2, label='AR(2)')
    plt.plot(smpl_range, list_rls, label='RLS(S-M)')
    plt.plot(smpl_range, list_gl, label='glasso')
    plt.legend()
    plt.xticks(smpl_range[::num_ranges//10])
    # plt.yscale('log')
    plt.grid()
    plt.show()
    
def generate_normal(num_samples):
    return np.random.normal(0, 1, num_samples)

def generate_corr_normal(num_samples):
    return normalize(make_1st_correlation(generate_normal(num_samples), 0.95))
In [19]:
# i.i.d.雑音に対して試す
plot_inv_auto_corr_rmse_curve(4, 0.999, np.arange(50, 1000, 20), generate_normal, 20)
plot_inv_auto_corr_rmse_curve(16, 0.999, np.arange(50, 1000, 20), generate_normal, 10)
In [20]:
# 相関付き雑音に対して試す
plot_inv_auto_corr_rmse_curve(4, 0.999, np.arange(50, 1000, 20), generate_corr_normal, 20)
plot_inv_auto_corr_rmse_curve(16, 0.999, np.arange(50, 1000, 20), generate_corr_normal, 10)

考察

  • AR(1), AR(2)は大数の法則的に漸近しているように見える（1/√n）。
    • これは最尤推定量を使ったからと想像。指数移動平均でやっても結果はほぼ同じだった。
    • これだと遅い、ということなのかもしれない。
      • 最初の推定値が悪すぎるらしい。それによってRLSと差が出てしまっているっぽい。
      • RLSの推定値が "良い" とはどういうことなのだろうか。。。
    • AR(2)の方がわずかに残差が小さい傾向。
  • glassoはARよりちょっと良いくらいか？相関付き雑音だと却って悪いくらいの性能。
    • 正則化パラメータを小さくすると良くなるが、相関付き雑音で条件数が悪くなって警告・エラーを出し始める

次確かめること

  • 大域的な真値への収束特性確認。AR(1)が良かった記憶があるが、果たして本当にそうだったか。
In [21]:
# 真の逆行列（固定）との誤差RMSEの変化をプロット
def plot_fixed_inv_auto_corr_rmse_curve(dim, forget, smpl_range, data_generator, num_trials):
    num_ranges = smpl_range.size
    list_ar1 = np.zeros(num_ranges)
    list_ar2 = np.zeros(num_ranges)
    list_rls = np.zeros(num_ranges)
    list_gl = np.zeros(num_ranges)
    for _ in range(num_trials):
        # 長期サンプルにおける真の逆行列を先に求めておく
        wholedata = data_generator(np.max(smpl_range))
        S = calc_auto_corration_matrix(wholedata, dim)
        Sinv = np.linalg.inv(S)
        Sinv /= np.linalg.norm(Sinv)
        for idx, nsmpl in enumerate(smpl_range):
            data = wholedata[0:nsmpl]
            Sinv_ar1 = calc_inv_auto_corr_by_ar1(data, dim, forget)
            Sinv_ar2 = calc_inv_auto_corr_by_ar2(data, dim, forget)
            _, Sinv_rls = calc_inv_auto_corr_by_rls(data, dim, forget)
            # Sinv_gl = calc_inv_auto_corr_by_gl(data, dim, 1e-1)
            Sinv_ar1 /= np.linalg.norm(Sinv_ar1)
            Sinv_ar2 /= np.linalg.norm(Sinv_ar2)
            Sinv_rls /= np.linalg.norm(Sinv_rls)
            # Sinv_gl /= np.linalg.norm(Sinv_gl)
            list_ar1[idx] += np.linalg.norm(Sinv - Sinv_ar1)
            list_ar2[idx] += np.linalg.norm(Sinv - Sinv_ar2)
            list_rls[idx] += np.linalg.norm(Sinv - Sinv_rls)
            # list_gl[idx] += np.linalg.norm(Sinv - Sinv_gl)
    list_ar1 /= num_trials
    list_ar2 /= num_trials
    list_rls /= num_trials
    # list_gl /= num_trials
    plt.plot(smpl_range, list_ar1, label='AR(1)')
    plt.plot(smpl_range, list_ar2, label='AR(2)')
    plt.plot(smpl_range, list_rls, label='RLS(S-M)')
    # plt.plot(smpl_range, list_gl, label='glasso')
    plt.legend()
    plt.xticks(smpl_range[::num_ranges//10])
    # plt.yscale('log')
    plt.grid()
    plt.show()
In [22]:
# i.i.d.雑音に対して試す
plot_fixed_inv_auto_corr_rmse_curve(4, 0.999, np.arange(20, 700, 10), generate_normal, 20)
plot_fixed_inv_auto_corr_rmse_curve(8, 0.999, np.arange(20, 700, 10), generate_normal, 20)
plot_fixed_inv_auto_corr_rmse_curve(16, 0.999, np.arange(20, 700, 10), generate_normal, 20)
In [23]:
# 相関付き雑音に対して試す
plot_fixed_inv_auto_corr_rmse_curve(4, 0.999, np.arange(20, 700, 10), generate_corr_normal, 20)
plot_fixed_inv_auto_corr_rmse_curve(8, 0.999, np.arange(20, 700, 10), generate_corr_normal, 20)
plot_fixed_inv_auto_corr_rmse_curve(16, 0.999, np.arange(20, 700, 10), generate_corr_normal, 20)

考察

  • 開始直後の収束は早い。途中でRLS(Sherman-Morrison)に追い抜かれる
  • 次数が大きいと定常状態にすぐ入ってしまう。そしてS-Mに追い抜かれる。
    • でも次数4ならばAR(2)は喰らいついていく感じがある。結局行列要素全部求めてるから当然か。
  • AR(1), AR(2)仮定、結局よくなさそう、というのでクローズしそう。

次は

  • 共役勾配法を実装し直す

報告用プロット

In [24]:
plot_auto_corr_compare_ar1_and_gt(gaussnoise, 16, 0.999)
plt.savefig("InverseAutoCorrCompareAR1_iid.png")
plot_auto_corr_compare_ar1_and_gt(corr_gaussnoise, 16, 0.999)
plt.savefig("InverseAutoCorrCompareAR1_corr.png")

追加実験: 実データの自己相関行列の逆

In [54]:
import soundfile

realdata, _ = soundfile.read('./data/seijya_no_kousin_head4s_Lch.wav')

for smpl in np.arange(5000, 50000, 5000):
    partial_data = realdata[smpl-5000:smpl]
    print('%d-%d samples:' % (smpl-5000, smpl))
    plot_auto_corr_compare_ar1_and_gt(partial_data, 16, 0.999)
    plot_auto_corr_compare_ar2_and_gt(partial_data, 16, 0.999)
    plot_auto_corr_compare_gl_and_gt(partial_data, 16, 0.01)
    plt.show()

0-5000 samples:

5000-10000 samples:

10000-15000 samples:

15000-20000 samples:

20000-25000 samples:

25000-30000 samples:

30000-35000 samples:

35000-40000 samples:

40000-45000 samples:

考察

  • 実データの自己相関行列がほぼ変わらない
    • 真値:行列の中央が最大値(約0.2)で、端に行くほど減衰。
    • 真値:微妙に波打ってる。
  • 真値、AR(1)はずっと変わってない
    • (1,1), (N,N)で大きな差異が出ているのが気になる。
  • glassoは僅かに変わっているが、対角優位で端に行くほど減衰する傾向は変わらない
    • AR(1):glassoの推定結果に近い
In [ ]: